####Permutation tests####
#Remove objects
rm(list = ls())

#Load pkgs
library(tidyverse)
library(parallel)

#Load helper functions
source("scripts/helper_functions.R")

#Load main dataset
df_merged_mean <- read_delim("input/df_merged_mean.csv", 
                             ";", escape_double = FALSE, col_types = cols(date_text = col_date(format = "%d/%m/%Y")), 
                             trim_ws = TRUE)

#Read positively related topics from main model and category names (necessary to run analysis.R before)

pos_names <- readRDS("output/pos_names.rds")
names_output <- readRDS("output/names_output.rds")
models_4 <- readRDS("output/model_4.rds")

#Permutation Tests
df_merged_mean$X1 <- seq(3743)

df_merged_mean <- df_merged_mean %>% ungroup()

perm_test <- function(j){
  library(tidyverse)
  library(brglm)
  myTryCatch <- function(expr) {
    warn <- err <- NULL
    value <- withCallingHandlers(
      tryCatch(expr, error=function(e) {
        err <<- e
        NULL
      }), warning=function(w) {
        warn <<- w
        invokeRestart("muffleWarning")
      })
    list(value=value, warning=warn, error=err)
  }
  
for (i in seq(1000)){ #Number of runs changed from 10000 (as in manuscript) to 1000 to run faster
  set.seed(123+i) #seeds changes depending on run number to generate different combinations
  print(i)
  df_merged_mean_temp <-df_merged_mean %>% mutate(
   off_rob = base::sample(off_rob))
  df_merged_mean_temp <- df_merged_mean_temp %>%  group_by(date_text) %>%
      mutate(var_time_date= var(off_rob)) %>%
      mutate(date_pos= ifelse(var_time_date>0,date_text,0)) %>% ungroup()
  df_merged_mean_temp <- df_merged_mean_temp %>%
  group_by(website_id) %>%
    mutate(off_rob_l= dplyr::lag(off_rob,order_by=website_id)) %>% ungroup()
  e1 <- myTryCatch(model <- brglm(off_rob~scale(get(j))+off_rob_l+factor(website_id)+factor(date_pos),data=df_merged_mean_temp,
               family=binomial(link="logit"),method="brglm.fit",pl=T))
  if(e1$warning[1]=="Iteration limit reached"){ 
  coeff <- "NA"
    if (i==1){
      coeff1 <- coeff
    }
    if (i>1){
      coeff1 <- rbind(coeff1,coeff)
    }}
  else{
    coeff <- model$coefficients[2]
    if (i==1){
      coeff1 <- coeff
    }
    if (i>1){
      coeff1 <- rbind(coeff1,coeff)
    }
  }
} #Randomly shuffles 19 attacks days around and runs models
return(coeff1)
}

run_perm_test <- function(j){
  # Run marginal models
  start.time <- Sys.time()
  library(parallel)
  
  # Calculate the number of cores
  no_cores <- detectCores(all.tests = FALSE, logical = TRUE)
  no_cores <- no_cores-1
  
  # Initiate cluster
  cl <- makeCluster(no_cores)
  
  clusterExport(cl, list("df_merged_mean","perm_test","myTryCatch"))
  
  marginal_models <- parLapply(cl, 
                               c(j), 
                               perm_test)
  
  stopCluster(cl)
  end.time <- Sys.time()
  return(marginal_models)
}

#Simulated estimates (runs approx. 6h when simulations are set to 10,000)
perm_sims <- run_perm_test(pos_names)

perm_sims_0 <- perm_sims

#Observed estimates
for (i in pos_names){
  pos_names_estimates <- models_4[1,][[i]]$coefficients[2]
  if (i == "opposition_exiled"){
    pos_names_out <- pos_names_estimates  
  }
  else
    pos_names_out <- c(pos_names_out,pos_names_estimates)
}


pdf("output/figures/Figure_D1_permutation.pdf",width=16, height=8)
par(mfrow=c(3,4))
for (i in seq(1:11)){
  plot(density(as.numeric(perm_sims[[i]]),na.rm=T),main=
         as.character(filter(names_output,names_fe %in% pos_names)$topic_names_text_fe)[i],
       xlim = c(-1.5, 1.5))
  abline(v=models_4[1,][[pos_names[i]]]$coefficients[2],col="red")
  abline(v=quantile(as.numeric(perm_sims[[i]]),0.95,na.rm=T),lty=2)
}
dev.off()


